#include <stdlib.h>
#include <StandardStream.h>
#include <OPS_Globals.h>
#include <ArrayOfTaggedObjects.h>
#include <OPS_Stream.h>
#include <DOF_Numberer.h>
#include <Domain.h>
#include <ID.h>
#include <Node.h>
#include <LinearCrdTransf3d.h>
#include <PDeltaCrdTransf3d.h>
#include <CorotCrdTransf3d.h>
#include <ElasticMaterial.h>
#include <ElasticPPMaterial.h>
#include <Steel01.h>
#include <ElasticIsotropicThreeDimensional.h>
#include <DruckerPrager3D.h>
#include <PlateFiberMaterial.h>
#include <Fiber.h>
#include <UniaxialFiber3d.h>
#include <FiberSection3d.h>
#include <ElasticSection3d.h>
#include <SectionAggregator.h>
#include <LobattoBeamIntegration.h>
#include <FixedLocationBeamIntegration.h>
#include <ElasticBeam3d.h>

StandardStream sserr;
OPS_Stream *opserrPtr = &sserr;
int main(int argc, char **argv) {

Domain *theDomain = new Domain();
int i = 1;
int j = 1;

Node *node1 = new Node(1, 6, 0.000000, 0.000000, 0.000000);
theDomain->addNode(node1);
Node *node2 = new Node(2, 6, 1.000000, 0.000000, 0.000000);
theDomain->addNode(node2);
Node *node3 = new Node(3, 6, 0.000000, 0.000000, 2.000000);
theDomain->addNode(node3);

Vector Trans0(3);
Trans0(0) = 0.0;
Trans0(1) = 0.0;
Trans0(2) = 0.0;
Vector XZ1(3);
XZ1(0) = 1.000000;
XZ1(1) = 0.000000;
XZ1(2) = 0.000000;
CrdTransf *Corotational1 = new CorotCrdTransf3d(1, XZ1, Trans0, Trans0);
Vector XZ2(3);
XZ2(0) = 1.000000;
XZ2(1) = 0.000000;
XZ2(2) = 0.000000;
CrdTransf *PDelta1 = new PDeltaCrdTransf3d(2, XZ2, Trans0, Trans0);

UniaxialMaterial *ElasticMaterial1 = new ElasticMaterial(1, 205000.000000);
UniaxialMaterial *ElasticMaterial2 = new ElasticMaterial(2, 206000.000000);
UniaxialMaterial *ElasticPPMaterial1 = new ElasticPPMaterial(3, 205000.000000, 0.001000);
UniaxialMaterial *ElasticPPMaterial2 = new ElasticPPMaterial(4, 205000.000000, 0.001000, -0.000500, 0.000000);
UniaxialMaterial *ElasticPPMaterial3 = new ElasticPPMaterial(5, 205000.000000, 0.001000, -0.001000, 0.001000);
UniaxialMaterial *Steel01Material1 = new Steel01(6, 345.000000, 205000.000000, 0.010000, 0.000000, 1.000000, 0.000000, 1.000000);

NDMaterial *ElasticIsotropic1 = new ElasticIsotropicThreeDimensional(7, 205000.000000, 0.200000, 0.000000);
NDMaterial *DruckerPrager1 = new DruckerPrager3D(8, 185000.000000, 90000.000000, 345.000000, 0.000000, 0.000000, 1.000000, 1.000000, 0.100000, 0.100000, 1850.000000, 0.100000, 0.000000, 101.000000);
NDMaterial *DruckerPrager2 = new DruckerPrager3D(9, 185000.000000, 90000.000000, 390.000000, 0.000000, 0.000000, 1.000000, 1.000000, 0.100000, 0.100000, 1850.000000, 0.200000, 0.000000, 101.000000);
NDMaterial *ElasticIsotropic1Plane = new PlateFiberMaterial(1007, *ElasticIsotropic1);
NDMaterial *DruckerPrager1Plane = new PlateFiberMaterial(1008, *DruckerPrager1);
NDMaterial *DruckerPrager2Plane = new PlateFiberMaterial(1009, *DruckerPrager2);

double locy;
double locz;
double AFiber;
int tempno;
Fiber *fibers1[180];
tempno = 0;
for (i = 1; i <= 10; i++) {
   for (j = 1; j <= 9; j++) {
      locy = -10.000000 + 20.000000 / 10 * 0.5 + (i - 1) * 20.000000 / 10;
      locz = -20.000000 + 40.000000 / 9 * 0.5 + (j - 1) * 40.000000 / 9;
      AFiber = 20.000000 * 40.000000 / (10 * 9);
      fibers1[tempno] = new UniaxialFiber3d(tempno + 1, *ElasticPPMaterial1, AFiber, locy, locz);
      tempno = tempno + 1;
   }
}
for (i = 1; i <= 10; i++) {
   for (j = 1; j <= 9; j++) {
      locy = -30.000000 + 40.000000 / 10 * 0.5 + (i - 1) * 40.000000 / 10;
      locz = -40.000000 + 60.000000 / 9 * 0.5 + (j - 1) * 60.000000 / 9;
      AFiber = 40.000000 * 60.000000 / (10 * 9);
      fibers1[tempno] = new UniaxialFiber3d(tempno + 1, *Steel01Material1, AFiber, locy, locz);
      tempno = tempno + 1;
   }
}
Fiber *fibers2[115];
tempno = 0;
for (i = 1; i <= 10; i++) {
   for (j = 1; j <= 9; j++) {
      locy = -30.000000 + 40.000000 / 10 * 0.5 + (i - 1) * 40.000000 / 10;
      locz = -40.000000 + 60.000000 / 9 * 0.5 + (j - 1) * 60.000000 / 9;
      AFiber = 40.000000 * 60.000000 / (10 * 9);
      fibers2[tempno] = new UniaxialFiber3d(tempno + 1, *Steel01Material1, AFiber, locy, locz);
      tempno = tempno + 1;
   }
}
for (i = 1; i <= 5; i++) {
   for (j = 1; j <= 5; j++) {
      locy = 10.000000 + 20.000000 / 5 * 0.5 + (i - 1) * 20.000000 / 5;
      locz = 20.000000 + 20.000000 / 5 * 0.5 + (j - 1) * 20.000000 / 5;
      AFiber = 20.000000 * 20.000000 / (5 * 5);
      fibers2[tempno] = new UniaxialFiber3d(tempno + 1, *ElasticPPMaterial3, AFiber, locy, locz);
      tempno = tempno + 1;
   }
}
SectionForceDeformation *FiberSection1[10];
FiberSection1[0] = new FiberSection3d(1, 180, fibers1);
FiberSection1[1] = new FiberSection3d(2, 180, fibers1);
FiberSection1[2] = new FiberSection3d(3, 180, fibers1);
FiberSection1[3] = new FiberSection3d(4, 180, fibers1);
FiberSection1[4] = new FiberSection3d(5, 180, fibers1);
FiberSection1[5] = new FiberSection3d(6, 180, fibers1);
FiberSection1[6] = new FiberSection3d(7, 180, fibers1);
FiberSection1[7] = new FiberSection3d(8, 180, fibers1);
FiberSection1[8] = new FiberSection3d(9, 180, fibers1);
FiberSection1[9] = new FiberSection3d(10, 180, fibers1);
SectionForceDeformation *FiberSection2[7];
FiberSection2[0] = new FiberSection3d(1, 180, fibers1);
FiberSection2[1] = new FiberSection3d(2, 180, fibers1);
FiberSection2[2] = new FiberSection3d(3, 180, fibers1);
FiberSection2[3] = new FiberSection3d(4, 115, fibers2);
FiberSection2[4] = new FiberSection3d(5, 115, fibers2);
FiberSection2[5] = new FiberSection3d(6, 115, fibers2);
FiberSection2[6] = new FiberSection3d(7, 115, fibers2);
SectionForceDeformation *ElasticSection1[7];
for (i = 0; i < 7; i++)
ElasticSection1[i] = ElasticSection3d(3, 205000.000000, 0.100000, 0.010000, 0.010000, 105001.000000, 0.020000);
UniaxialMaterial *MatCombinedForSA1[3];
MatCombinedForSA1[0] = ElasticMaterial1;
MatCombinedForSA1[1] = ElasticPPMaterial2;
MatCombinedForSA1[2] = Steel01Material1;
ID DofsForSA1(3);
DofsForSA1(0) = 2;
DofsForSA1(1) = 6;
DofsForSA1(2) = 4;
SectionForceDeformation *AggregatedSection1 = new SectionAggregator(4, *FiberSection2[0], 3, &MatCombinedForSA1[0], DofsForSA1);
SectionForceDeformation *AggregatedSections1[7];
AggregatedSections1[0] = AggregatedSection1;
AggregatedSections1[1] = AggregatedSection1;
AggregatedSections1[2] = FiberSection2[0];
AggregatedSections1[3] = FiberSection2[0];
AggregatedSections1[4] = FiberSection2[0];
AggregatedSections1[5] = AggregatedSection1;
AggregatedSections1[6] = AggregatedSection1;

BeamIntegration *LobattoIntegration1 = new LobattoBeamIntegration();
Vector VecLoc1(7);
for (i = 0; i < 7; i++)
VecLoc1(i) = 1.0/6*(double)i;
BeamIntegration *FixedLocationIntegration1 = new FixedLocationBeamIntegration(7, VecLoc1);
Vector VecLoc2(7);
VecLoc2(0) = 0.000000;
VecLoc2(1) = 0.100000;
VecLoc2(2) = 0.200000;
VecLoc2(3) = 0.300000;
VecLoc2(4) = 0.700000;
VecLoc2(5) = 0.800000;
VecLoc2(6) = 1.000000;
BeamIntegration *FixedLocationIntegration2 = new FixedLocationBeamIntegration(7, VecLoc2);

ElasticBeam3d *ElasticBeamColumn1 = new ElasticBeam3d(1, 0.100000, 205000.000000, 109000.000000, 0.020000, 0.010000, 0.010000, 1, 2, *Corotational1);
ElasticBeam3d *ElasticBeamColumn2 = new ElasticBeam3d(2, 0.100000, 205000.000000, 109000.000000, 0.020000, 0.010000, 0.010000, 2, 3, *PDelta1);


}